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-^The  scope  of  this  study  is  to  develop  the  basic  equations  for  deriving  a  finite  element  for¬ 
mulation  which  can  be  used  to  solve  problems  related  to  convection  and  diffusion  dominated  flows. 

The  formulation  is  based  on  the  introduction  of  a  generalized  quantity  defined  as  heat 
displacement.  The  governing  equation  is  expressed  in  terms  of  this  quantity  and  a  variational 
formulation  is  developed  which  leads  to  an  equation  similar  in  form  to  Lagrange's  equation  of 
mechanics.  This  equation  may  be  solved  by  any  numerical  method,  though  it  is  of  particular  . _ 
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for  application  of  the  finite  element  method. 

The  developed  formulation  is  used  to  derive  two  finite  element  models  for  solving  convection- 
diffusion  boundary  value  problems.  The  performance  of  the  two  element  models  is  investigated  and 
numerical  results  are  given  for  different  cases  of  convection  and  diffusion  with  three  types  of  boundary 
conditions. 

The  numerical  results  obtained  show  not  only  the  efficiency  of  the  numerical  models  to  handle 
pure  convection,  pure  diffusion  and  mixed  convection-diffusion  problems  but  also  good  stability  and 
accuracy.  The  applications  of  the  developed  numerical  models  is  not  limited  to  diffusion-convection 
problems  but  can  also  be  applied  to  other  types  of  problems  such  as  mass  transfer,  hydrodynamics  and 
wave  propagation.  ,. 
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FINITE  ELEMENT  MODELING  FOR 
CONVECTION-DIFFUSION  PROBLEMS 


INTRODUCTION 

The  simulation  of  thermally  induced  waves  requires  the  solution  of  the  convection- 
diffusion  equation.  Analytical  and  numerical  solutions  of  this  equation  have  attracted  consider¬ 
able  attention  in  a  variety  of  engineering  fields  due  to  its  wide  applicability. 

The  theory  for  convective  or  diffusive  dominated  flows  has  been  well-established  and  a 
variety  of  classical  approaches  exist  in  the  literature  for  the  solution  of  problems  in  this  area. 
One  such  solution  is  given  by  Price,  Cavendish  and  Varga  [1],  Analytical  solutions  are  valid 
primarily  for  linear  equations  and  their  application  to  problems  of  practical  interest  presents 
difficulties  due  to  the  limitations  of  such  solutions.  It  is  because  of  the  restrictive  nature  of  the 
analytical  solutions  that  research  efforts  have  been  focused  on  approximate  or  numerical  solu¬ 
tions  of  the  convection-diffusion  equation.  A  review  and  comparison  of  available  numerical 
methods  can  be  found  in  Lee  et  al  (1976),  Ehlig  (1977),  and  Genuchten  (1977).  Numerical 
methods  discussed  in  these  papers  include  finite  differences  and  some  finite  element  approxima¬ 
tions  15].  Most  of  the  numerical  schemes  produce  results  of  acceptable  accuracy  either  for  con¬ 
vection  dominated  flows,  or  for  diffusion  dominated  flows,  but  they  lack  uniformity  in  perfor¬ 
mance.  Finite  difference  schemes  are  the  least  attractive  ones  due  to  their  instability,  large 
oscillations,  and,  for  some  of  them  inherent  artificial  diffusion.  On  the  other  hand  finite  ele¬ 
ment  schemes  have  produced  more  reliable  numerical  solutions  but  their  application  is  limited 
to  certain  types  of  dispersion  while  performing  poorly  for  other  types.  Another  disadvantage  of 
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existing  finite  element  solutions  is  that  they  are  based  on  formulations  restricted  by  the  condi¬ 
tions  of  a  particular  problem  and  they  are  limited  in  applications  to  other  types  of  problems. 

In  view  of  the  limitations  and  lack  of  uniformity  of  existing  approaches,  it  is  desirable  to 
develop  a  unified  formulation  which  is  capable  of  treating  not  only  pure  convection-through- 
dispersion  to  pure  diffusion  problems,  but  also  other  types  of  problems  governed  by  similar 
equations.  In  order  to  achieve  this,  a  unified  variational  formulation  is  introduced  for  the 
convection-diffusion  equation;  this  equation  is  written  in  terms  of  a  generalized  quantity, 
defined  as  heat  displacement  (6,7],  which  is  similar  to  a  mechanical  displacement  and  has  units 
of  length.  As  of  this  definition,  changes  in  temperature  are  treated  as  thermal  deformations 
which  are  similar  to  mechanical  strains. 

In  the  first  part  of  this  study,  the  basic  definitions  are  introduced  and  the  convection- 
diffusion  equation  is  expressed  in  terms  of  the  heat  displacement.  A  variational  formulation  is 
then  derived,  based  on  the  principle  of  virtual  work  in  mechanics,  and  by  using  generalized 
coordinates  the  variational  equation  is  written  in  a  form  equivalent  to  that  of  the  Lagrangian 
equation  in  mechanics.  Since  the  derived  equation  is  expressed  in  terms  of  generalized  coordi¬ 
nates,  it  is  applicable  to  a  wide  variety  of  physical  problems  and  can  be  solved  by  any  numerical 
method. 

The  variational  form  of  the  derived  equation  is  most  suitable  for  applying  the  finite  ele¬ 
ment  method  for  its  numerical  solution.  This  is  done  in  the  second  part  of  this  study,  where 
the  basic  finite  element  method  is  used  to  derive  two  finite  element  models  for  solving  initial  or 
boundary  value  problems.  The  first  model  is  based  on  a  linear  approximation  of  the  displace¬ 
ment  and  the  second  on  a  third  order  approximation.  The  matrix  equation  for  the  linear  model 
is  expressed  in  terms  of  nodal  displacements  and  for  the  higher  order  model,  in  terms  of  nodal 
displacements  and  nodal  temperatures. 
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Each  element  model  can  be  used  for:  a)  pure  convection,  b)  pure  diffusion,  and  c)  mixed 
convection-diffusion  problems  with  the  appropriate  boundary  conditions.  Also,  dynamic  or 
quasi-static  solutions  can  be  obtained  respectively  by  retaining  or  neglecting  the  inertia  term. 

Numerical  results  are  given  in  the  third  part  of  this  study  where  a  third-order,  backward 
finite  difference  scheme  is  employed  for  the  solution  of  the  system  of  differential  equations 
18,9].  For  the  combined  finite  element-finite  difference  scheme  the  stability,  convergence,  and 
accuracy  are  investigated  and  some  uniform  convergence  criteria  are  discussed.  Numerical 
results  are  also  given  for  a  number  of  convection  through  diffusion  cases  and  for  three  types  of 
boundary  conditions.  The  present  results  are  compared  to  existing  analyical  solutions  and  the 
accuracies  of  the  two  finite  element  models  are  discussed. 


1.  BASIC  EQUATIONS 

Consider  an  incompressible  medium  in  a  flow  field  subjected  to  external  heating.  Intially 
the  medium  is  at  a  uniform  temperature  T0,  which  will  be  referred  to  as  the  reference  tempera¬ 
ture,  and  the  state  at  this  temperature  will  be  referred  to  as  the  reference  state. 


The  instantaneous  absolute  temperature  is  denoted  by  T,  and  the  difference  T—  T0  defines 
the  instantaneous  relative  temperature  A0,  which  is  a  function  of  the  space  coordinates  and 
time.  Let 


T-T. 


M. 

T 

1  ft 


(1) 


be  defined  as  the  temperature  change  per  unit  temperatue  T0,  or  the  instantaneous  relative 
temperature  per  unit  temperature.  In  the  following  it  will  be  referred  to  as  the  temperature  9 
or  the  dimensionless  temperature. 


Assuming  cartesian  coordinates  ( x ,,  i  -  1,2,3),  the  temperature  field  9  satisfies  the 


convection-diffusion  equation 
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«f&!l  +  J_  ik,»  bl,» 

3r  3-t,  3x,  y  3.*} 

where  V,  is  the  velocity  vector  and  ku  the  thermal  diffusivity  (ktJ  -* 


We  now  define  a  vector  field  //,  ( xJti )  as  the  heat  displacement  vector  such  that 


The  summation  convention  is  assumed  for  repeated  indices  throughout  this  study. 


In  the  above  definition,  Eq.  (3),  0  represents  a  thermal  strain  analogous  to  mechanical 
strain.  Note  that  H,(xj,t)  has  the  dimension  of  displacement,  which  makes  it  analogous  to  a 
mechanical  displacement.  Thus  there  is  a  one-to-one  correspondence  between  heat 
displacement-mechanical  displacement  and  temperature  strain. 


Using  the  definition  from  Eq.  (3),  we  now  write  Eq.  (2)  as  follows 


Ml  +ve-k  - 

3 1  'e  k,J  dXj 


+  k,j  Vk 


cT0  3  Xj 


where  the  thermal  stress  <r  is  defined  by 


<r  •  c  T0  9  (6) 

with  c  the  heat  capacity  per  unit  volume,  and  k,j  »  (ku)~]. 

In  the  above  analysis  the  thermal  flow  field  is  governed  by  the  three  Eqs.  (3),  (5)  and  (6) 
which,  together  with  the  appropriate  boundary  conditions,  provide  a  complete  formulation  for 
convective  heat  transfer.  They  are  analogous  to  the  kinematic  relations,  stress-strain  relations 
and  momentum  equations  in  mechanics. 

The  advantage  of  introducing  the  heat  displacement  vector  are  more  apparent  when  the 
concept  of  virtual  work  is  used  to  derive  the  variational  formulation.  Furthermore,  the  above 
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analysis  is  more  suitable  for  applying  the  finite  element  method  which  is  based  on  displacement 
approximations. 


2.  VARIATIONAL  FORMULATION 


Following  the  usual  procedures  of  the  principle  of  virtual  work  in  mechanics,  we  consider 
that  the  medium  is  subjected  to  an  arbitrary  infinitesimal  virtual  displacement  8//,  from  the 
equilibrium  configuration.  The  corresponding  variations  80  are  given  by  £q.  (3).  Multiplica¬ 
tion  of  Eq.  (4)  by  f>H,  and  integation  over  a  volume  v  of  the  medium  yields 


J. 


—■  +  V,0  -  k  f>H,dv  «  0 

bt  J  bx, 


Integrating  by  parts  and  applying  the  divergence  theorem,  one  obtains 


X  *H,dv  +  X  mu.dv  +  X  M  <*/*,><*  =  X  W*HlVjdS 

where  n ;  is  the  unit  normal  vector  pointing  outward  at  the  boundary  surface  S.  From  Eq.  (3) 
one  derives 

XM  —  (hH,)dv  -  Jbfjk^Odv  =  k„f>E  (7) 

where  8y  is  the  Kronecker’s  delta  and  the  scalar  E  is  defined  as 

£  -  4-  f  02  dv  (8) 

2  v 

and  plays  the  role  of  a  potential  function.  Equation  (6)  is  now  written  as  follows 


/d  H  t*  c 

~~  8 Hjdv  +  Jk, j  VfihHjdv  -  Xs  KHfljdS.  (9) 

Eq.  (9)  may  be  considered  as  a  variational  principle  in  a  broad  sense  and  a  more  compact  form 

of  this  equation  can  be  derived  by  introducing  generalized  coordinates  defined  as 


Hj(Xj,t)  -  Hj(qv,Xj,t)  (10) 

where  the  generalized  coordinates  q„  are  functions  of  time.  The  advantage  of  using  generalized 

coordinated  is  that  H{  may  be  expressed  in  different  functional  forms.  Care  should  be  taken 

when  the  time  derivative  of  H,  is  considered,  and  it  should  be  expressed  as 
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H.  = 


dH, 


Q,  + 


a/^ 

dt  ' 


(11) 


In  terms  of  arbitary  variations  of  the  generalized  coordinates,  the  corresponding  displace¬ 
ment  field  variations  are 


8  //, 

f>H,=  — —  S  q, 


and  the  variation  of  the  potential  E  becomes 


S£  =  hq, 

H 


(12) 


(13) 


In  view  of  Eqs.  (12)  and  (13),  Eq.  (9)  may  be  written  for  each  arbitrary  variation  ft./*  as  follows 

|£  +  f  *  Ml  M.  *  +  /  x  Vs.  ,14) 

8  a,  Jv  IJ  dt  da,  1  da,  dqk  J 


dqk  IJ  dt  dqk 

From  Eq.  (11)  one  derives 


dH,  dH, 
d  qk  dqk 

and  the  second  term  in  Eq.  (14)  is  then  expressed 

dD 
dqk 

where 


f  dH,  d Hj  r  dH,  dH  0 

X^ir  "  XXtf  it  v*  -  ^ 


tX  KjH,Hjdy\ 


D  =  {X  X*  dv 

Eq.  (14)  then  takes  the  form 


where 


dD  ,  dE 
dqk  9qk 


+ 


Lk 


Qk 


and 


dH. 

K 'ft  ~  dv 


(15) 


(16) 


(17) 


Qk  -  f#v, 


dH, 

dqk 


dS 
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The  quaniily  ft  can  be  regarded  as  the  thermal  force  applied  on  the  boundary  S  of  the 
medium  of  volume  v.  The  quantity  Lk  represents  the  convection  term  and  for  the  case  of  a 
solid  medium  ( V,  =0)  it  is  equal  zero  and  Eq.  (16)  reduces  to  a  form  similar  to  Lagrange's 
equation  in  mechanics.  Eq.  (16)  as  it  was  derived  is  quite  general  in  the  sense  that  it  can  be 
applied  to  different  types  of  media  with  different  material  properties. 


Consider  now  a  special  case  where  the  displacement  field  can  be  approximated  by  a  linear 
combination  of  the  generalized  coordinates  as  follows 

H,(Xj,t)  =  qk(t)fkl  (Xj)  k  =  1,  n,  j  —  1, 2, 3.  (18) 

In  Eq.  (18)  the  coefficients  qk(t)  represent  a  degree  of  freedom  and  the  function/*, (xt) 
specifies  the  extent  to  which  qk(i)  participates  in  the  function  H,(Xj,t).  In  finite  element 
analysis,  for  example,  Eq.  (18)  may  be  considered  as  the  distribution  function  of  the  displace¬ 
ment  field,  where  qk  can  then  be  taken  as  nodal  displacement  or  nodal  deformations  depending 
on  the  type  of  element  selected. 


Differentiating  Eq.  (18)  with  respect  to  time  and  space  we  obtain 


H,  -  qjk, 


(19) 


BH, 

9  ~  TT  Qkfku- 

Ox, 

The  scalar  £,  the  vector  Lk  and  the  invariant  D  from  Eqs.  (8),  (15),  and  (17)  are  expressed  in 
terms  of  Eqs.  (19)  as  follows 

E  *  2  ^mnQmQn 


^  “  2  dmnQmQn 


A k  “  8km  Qm 


(20) 
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where 

*mn  “  f 'mi.if nj.jdv 

dmn  *  ^  tjfmtfnj  dv 

«km  “  X  Xy  V'fkjfn,n.ndv  <21  ) 

Substituting  the  specific  forms  of  £,  A  and  Z.A  from  Eqs.  (20)  into  Eq.  (16)  one  obtains 

d.jQ,  +  ( qu  +  e,j)qj  =  0,  (22) 

with 

ft  =  X  «»,/„  ^  (23) 

Equations  (22)  constitute  a  system  of  /?  ordinary  differential  equations  for  the  unknown  field 

parameters  <?A  ( A  =  l,w),  and  it  may  represent  the  heat  displacement  field  H,.  This  system  of  n 
equations  can  be  solved  together  with  the  appropriate  boundary  condition  by  any  numerical 
technique.  Thus  this  variational  formulation  is  not  restricted  to  applications  of  the  finite  ele¬ 
ment  method  but  is  appropriate  for  applying  other  numerical  schemes  as  well. 

Another  advantage  of  the  derived  equations  is  that  they  are  not  restricted  to  solving 
convection-diffusion  problems.  By  appropriate  choice  of  the  variables  qk  to  represent  other 
physical  quantities,  the  derived  equations  can  be  used  to  solve  problems  involving  such  quanti¬ 
ties  as  concentration  or  velocity  fields. 

3.  FINITE  ELEMENT  ANALYSIS 

In  order  to  demonstrate  the  application  of  the  finite  element  method  to  the  previously 
derived  variational  formulation,  two  element  models  are  chosen  to  approximate  the  heat  dis¬ 
placement.  The  first  model  is  a  linear  element  with  minimum  degrees  of  freedom  (££)  and 
the  second  is  a  higher  order  element  with  four  degrees  of  freedom  (C£),  known  as  first  order 
cubic  Hermitian.  Although  both  elements  are  one  dimensional  approximations,  they  provide  a 
good  lest  case  for  the  performance  of  any  variational  formulations.  An  extension  into  the  two 
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dimensional  space  is  easily  obtained  since  the  derived  equations  are  of  general  form.  If  the 
heat  displacement  is  approximated  by  a  third  order  polynomial 

H(x,t)  =  an  +  fl).v  +  a2x2  +  a3x\ 
then  the  temperature  9  is  given  by 

Q(x,t)  =  a,  +  2a2x  +  3  a3x2 

where  a,  are  time  dependent  coefficients  to  be  determined  for  each  of  the  element  models 

a.  Linear  element  (LE) 

For  the  linear  element  of  length  /  the  conditions  at  the  nodal  points  are: 

at  x  =  0  —  H  =  H\ 
and 

at  x  =/—*//«  H2 

where  H\  and  H7  are  the  nodal  values  of  the  heat  displacement  and  the  coefficients  a,  are  given 
by 

a0=  //,(/),  a,  -  j(H2-  H,),  fl2=  o3«  0.  (26) 

Substituting  these  coefficients  in  Eqs.  (24)  and  (25)  yields 

//(x,f)=  1  -  yj  H]U)  +  ~H2(i)  (27) 

and 

9ix.t)  =  -  //,(/)). 

Note  that  within  each  element  9  varies  only  with  time  for  the  {LE)  approximations. 


(24) 

(25) 


The  matrix  coefficients  of  Eq.  (22)  are  evaluated  in  terms  of  Eqs.  (27)  as  follows 


8mn 


Av\- 1 
2k  “I 


(28) 
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and 


®  |,-o  “  ~  ~  A9-  Q,  y_t  =  0:  =  .40 

«h£rc  ,  is  ,h.  cross  sectional  area  of  the  demem.  *is  ,he  differs, iviiy  and  Frs  .he  fluid  »eio- 
city. 


In  terms  of  Eqs.  (28).  Eq.  (22)  yields 


6k 


1 

-I 


*2 


b.  Cubic  element  (CE) 


(29) 


I 


For  ,he  cubrc  Hermiiian  element  the  coefficiems  a,  are  evaluared  from  the  nodal  values  of 

H{* •"  and  SPa‘ial  deriva'ives  a*  U»  "«lef.  which  are  .he  nodal  values  of  .he  temperature 

O(x.t). 


The  conditions  at  the  nodes  are 


and 


where  «,.#.>  and  <*„».)  are  .he  nodal  values  of  the  he,,  displacement  and  iemperaiure 
respectively.  Thus,  the  coefficients  of  Eq.  (24)  can  be  found  as 

ft  i  *  Of  |  »  9 1 


a2-“  ^-((02  +  20,)/  -3(//2-  //,)] 

*  7  I(«2  +  *j)/  -  2{Hi  -  //,)] 
and  the  expressions  for  H(x,t)  and  00c./)  are  given  by 

H(x,t)  -  /,,</,  +  /l20j  +  +  /,4?4 

-  h |,9,  +  h\i<fi  +  +  A,494. 


(30) 


(31) 
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The  corresponding  e0,  dt],  and  gu  are  given  by 

36  3/  -36  3/ 

A  31  4 I2  -31  -l2 
e‘J  “  30/  -36  -31  36  -31 
31  -I2  -31  4I2 

156  22/  54  -13/ 

Alk  221  4/2  13/  — 3/2 

^*  420  54  13/  156  -22/ 

-13/  -3/2  -22/  4/2 

-210  42/  210  -42/ 

-42/  0  42/  — 7/2 

^  “  TlO  -210  -42/  210  42/ 

42/  7/2  -42/  0 

and  the  components  of  the  generalized  force  are 


Q\  “  —  A9 j,  O3  —  A9i,  Qj  —  04  —  0. 
In  terms  of  Eqs.  (27M30),  Eq.  (22)  yields 

"1  j«l  -0, 

W1  *  +  1W  +  Ml  Ji,  -  ® 

^2)  0 


(32) 


(33) 


(34) 


(35) 


(36) 


(37) 
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or 

\d)  |«|  +  IU)  +  Mil*)  «  l<?> 

For  the  solution  of  a  particular  problem,  the  finite  element  models  derived  above  are 
assembled  according  to  the  direct  stiffness  method  to  obtain  the  global  equations  The  formula¬ 
tion  of  the  overall  problem  is  not  complete  unless  boundary  conditions  are  taken  into  con¬ 
sideration.  The  system  of  n  equations  together  with  the  appropriate  boundary  conditions  can  be 
solved  by  any  numerical  technique  used  for  solving  ordinary  differential  equations.  In  the  fol¬ 
lowing  section,  the  above  system  of  equations  is  solved  for  two  types  of  boundary  conditions  by 
using  a  backward  differences  in  time  integration  technique. 

4.  BOUNDARY  VALUE  PROBLEM 

The  one  dimensional  case  of  the  convention-diffusion  equation  is  considered  here  to 
evaluate  the  two  finite  element  models  introduced  previously  with  the  following  initial  and 
boundary  conditions. 


a.  Convection  dominated  flow 

The  equation  to  be  solved  is 


dt  dx 


where  V  is  the  flow  velocity  of  constant  value. 
The  initial  conditions  are 


9(x. 0)  -  0,  0  <  *  L 
and  the  boundary  conditions  are 


I.  9 (O.r)  - 

II.  9  (0,r)  - 


LzLl 

T0 

TrTp 

To 

0 


t  >  0 

o  <  r  <  r„ 
to  <  t 


(38) 


(39) 


(40) 
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I.  Linear  element: 


\dc 1  l/fl  +  tl&J  +  le,,ll  [Hi  -  {0| 

with 


ki 


2  I 
1  2 


l&J-  3  K  K 


-1  1 
-1  1 


\et]  -  6  WjK0 


1  -I 
-1  1 


and  |0|  -  6 W0Ka 


-0 

9 


(45) 


II.  Cubic  element: 


with 


k) 


[H\ 


+  H&l  + 


(01 


165 

54 

22/ in 

-13/ 1^0 

54 

156 

13/ 

-22/ 

22/M^0 

13/  l*o 

4/  H'o 

-3/  ^ 

-13/  W0 

22/  W0 

-3/ 

4/^ 

kl 


2V0Wo 


-210 
—210 
— 42/  H'q 

42/ in 


210 
210 
42/  M^0 
-42 JW0 


42/  W0 
-42/  H^o 
0 

7 IWj 


-42/ in 
42/^0 
-7/lftf 
0 


(46) 


kJ-  l4WlK„ 


36 

-36 

3/ 

3/^0 

~0\ 

-36 

36 

-3/1^0 

-3/  in 

,  (0)  -  420  in*o 

3/^o 

-3/  W0  4 jWl 

-MWl 

0 

-3/H^o 

-3/^0  -1/^ 

4twl 

0 

where  W',  -  Z./4  and  and  AT,,  are  the  dimensionless  constant  velocity  and  diffusivity  respec¬ 
tively.  The  bar  over  the  variables  has  been  eliminated  for  simplicity.  The  boundary  conditions 
are  transformed  due  to  the  dimensionless  quantities  as  follows 
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a.  Convection  —dominated  flow 

1.0(0./)  -  1.0.  />0 

1 1.0.  0  <  /  ^  t0 

H.  0(0,/)  -  jo  t>to 

b.  Convection-diffusion 


1.0(0,/)  - 

1.0. 

/  > 

0 

00../)  - 

0. 

/  > 

0 

(l-o. 

0  < 

1  <  to 

II.  0(0./)  - 

|o.o. 

/  > 

to- 

0(L,t)  — 

0.0, 

/  ^ 

0. 

(47) 


(48) 


The  assembly  of  the  above  equations  for  the  overall  problem  and  their  modification  due 
to  the  boundary  conditions  is  coded  in  a  computer  program  given  in  Appendix  C.  After  solving 
for  the  displacements  of  the  (LE)  model,  the  temperture  for  the  /  '*  element  can  be  obtained 
through  the  relation 

e,-W0(Hl  +  l-Hl)i-\,n  (49) 

For  the  (CE)  model,  the  solution  of  the  system  of  equations  will  directly  give  nodal  displace¬ 
ments  as  well  as  nodal  temperatures. 


5.  NUMERICAL  SOLUTION 


The  one-dimensional  convection-diffusion  problem  has  been  formulated  by  the  finite  ele¬ 
ment  method  and  its  solution  can  be  obtained  from  the  system  of  ordinary  differential  equa¬ 
tions  in  matrix  form  presented  in  the  previous  section.  For  the  boundary  value  problem,  with 
given  boundary  conditions,  numerical  solutions  are  obtained  by  applying  suitable  numerical 
integration  techniques. 


IS 


ti  MHVVtltMS 


In  this  section,  the  numerical  errors  induced  hv  the  time  integration  scheme  and  by  the 
finite  element  method  are  evaluated  as  ate  the  convergences  of  the  solution  Furthermore, 
results  are  given  for  the  different  types  ol  boundary  conditions  considered  in  the  previous  sec¬ 
tion 


a.  Time  integration  technique 


The  systems  of  ordinary  differential  equations  obtained  for  the  two  element  models  are 
solved  by  using  a  third-order,  backward  finite  difference  approximation  The  general  formula 
for  the  first  derivative  can  be  written  in  the  following  form 


dy(t,)  |  jl 

-  4-  y  a.v  (/,_.)  +  0(A”  ')  (50) 

*  h  frQ  - 

where  h  is  the  size  of  the  lime  discretization  and  a ,  are  constant  coefficients  If  only  third  order 

and  lower  terms  are  retained  (n  **  4).  then  the  coefficients  are  given  by 


a„  - 


11.0 
60  ' 


a,  -  - 


18.0 

6.0  ' 


a :  — 


9.0 

6.0' 


o,  -  - 


2J0 

6.0' 


at  **  0.0 


(51) 


for  the  third-order  approximation. 


In  order  to  evaluate  the  stability  and  convergence  of  the  above  numerical  scheme,  results 
were  obtained  for  a  number  of  different  time-step  sizes.  A  stability  analysis  given  in  Appendix 
A  shows  that  the  scheme  is  unconditionally  stable.  This  stability  is  not  related  to  error  esti¬ 
mates  or  rate  of  convergence. 

Consider  the  case  of  pure  diffusion  (K0  —  1,  V0  -  0)  as  a  lest  case  for  illustrating  some 
numerical  results.  A  characteristic  length  L  —  5  is  chosen  to  represent  the  semi-infinite  space, 
which  is  divided  into  10  elements  (TNE  —  10).  Results  with  respect  to  time  are  given  for  the 
point  x  •  1.0  and  for  the  following  boundary  conditions 

0(0,/ > »  1.0,  r  £  0 

»(Lj)  -  0.0,  /  >  0 
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This  choice  of  boundary  conditions  is  appropriate  for  evaluating  numerical  solution  errors  and 
convergence.  Numerical  results  for  different  time  step  sizes  are  presented  in  Figs.  1,2  and 
Table  l  for  the  two  element  models.  The  optimum  time  step  was  found  to  be  within  the  range 
of  0.01  to  0.2.  For  smaller  or  larger  time  steps  the  error  of  the  numerical  solution  increases. 
Referring  to  Fig.  1,  it  is  apparent  that  as  the  time  step  becomes  smaller,  the  solution  does  not 
converge  monotonically.  Comparing  the  results  for  the  two  models,  the  (CE)  as  expected,  is 
much  more  accurate  than  the  (LE).  An  optimum  time  step  size,  approximately  in  the  range  of 
0.01  to  0.025,  yields  the  most  accurate  results. 

The  non-monotonic  convergence  is  due  to  the  fact  that  the  numerical  solution  depends 
not  only  on  the  time  step  size  but  also  on  the  spatial  discretization  as  well.  This  is  demon¬ 
strated  in  the  following,  where  numerical  experiments  provide  a  criterion  for  uniform  conver¬ 
gence. 

b.  Convergence  of  the  Finite  Element  Solution 

In  order  to  investigate  the  convergence  of  the  two  finite  element  models,  the  previously 
discussed  problem  is  considered  and  with  the  same  constants.  The  total  number  of  elements  by 
which  the  characteristic  length  L  is  represented  is  designated  as  TNE  and  the  numbers  of  ele¬ 
ments  from  x  —  0  to  x  —  1.0  is  denoted  by  NE.  The  results  of  this  part  are  given  for  two 
cases,  one  where  the  time  step  size  A  t  is  kept  constant  for  different  TNE,  and  the  other  case 
where  the  ratio  Ar/Ax2  is  kept  constant,  where  Ax  is  the  dimensionless  length  of  an  element. 

Convergence  For  Constant  A/ 

Results  for  this  case  are  given  for  the  time  step  size  A  t  —  0.025  for  both  element  models. 
For  the  (LE)  model  TNE  is  5,  10,  20,  30  with  the  NE  being  1,2,4  and  6  respectively.  For  the 
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Fig.  2— Temperature  distribution  at  x  —  1.0  as  a  function  of  time  for  the  (CE)  model,  TNE 
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Table  1.  Numerical  results  for  the  temperature  at  \  =  10  for 
the  (CE)  model  with  constant  TNE  =  10 


Time 

l> 

II 

© 

Nj 

© 

il 

A/  -  0  05 

A i  =  0  025 

Exact 

0.2 

0.050761 

0  061340 

0  0824^5  0  099961 

i 

0  113846 

0.4 

0.167021 

0  226432 

0  24h4.0 

0  254928 

0  263552 

0.6 

0.305074 

0  345081 

0  350841 

0  355827 

0  361310 

0.8 

0.405435 

0.414393 

0  421731 

0  425394 

0  429195 

1.0 

0.464985 

0.468150 

0  473867 

0.476695 

0  479500 

1.2 

0.503298 

0.509717 

0.514170 

0.516448 

0.518605 

1.4 

0.535127 

0.542831 

0.546491 

0.548389 

0.550090 

1.6 

0.563341 

0.570084 

0.573147 

0.57478 

0.576149 

1.8 

0.587542 

0.592999 

0.595609 

0.597047 

0.598161 

2.0 

0.608034 

0.612611 

0.614868 

0.616085 

0.617074 

2.2 

0.625647 

0.629639 

0.631614 

0.632687 

0.633553 

2.4 

0.641050 

0.644596 

0.646344 

0.647267 

0.648076 

(CE)  model  TNE  is  5,  10,  20,  25  with  NE  being  1,2, 4, 5  respectively.  In  Fig.  3  results  for  the 
(LE)  model  are  presented  for  the  temperature  at  x  =  1.0  with  respect  to  time  for  the  four 
different  values  of  TNE.  As  one  may  see  from  this  figure,  the  results  of  the  finite  element 
solution  do  not  converge  monotonically  to  the  exact  solution  as  TNE  increases. 

For  the  CE,  the  results  are  given  in  Table  2  for  the  time  history  of  the  temperature  of  x  — 
1.0  and  for  different  values  of  TNE  (i.e.  5,10,20,25).  By  increasing  the  TNE  the  results  show 
improvement.  However,  increasing  the  number  of  elements  does  not  necessarily  imply  uni¬ 
form  convergence. 

To  obtain  a  better  view  of  the  convergence  of  a  typical  data  point,  Fig.  4  shows  the  errors 
for  the  temperature  with  respect  to  TNE  for  constant  At  and  with  respect  to  At  for  constant 
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Fig.  3  — Temperature  distribution  at  x  =  1.0  as  a  function  of  time  for  the  (LE)  model,  A t  —  0.025. 
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Table  2.  Numerical  results  for  the  temperature  at  v  =  1.0  for 
the  (CE)  model  with  constant  At  =  0.025. 


Time 

TNE  -  5 

TNE  -  10 

TNE  -  20 

TNE  -  25 

Exact 

0.2 

0.072260 

0.099961 

0.102454 

0.102548 

|  0.113846 

0.4 

0.242260 

0.254928 

0.256015 

0.256053 

0.263552 

0.6 

0.350094 

0.355827 

0.356245 

0.356256 

0.361310 

0.8 

0.422881 

0.425394 

0.425537 

0.425540 

0.429195 

1.0 

0.475747 

0.476695 

| 

0.476705 

0.476708 

0.479500 

1.2 

0.516323 

0.516448 

0.516484 

0.516389 

0.518605 

1.4 

0.548726 

0.548389 

0.548282 

0.548282 

0.55090 

i 

1.6 

0.575355 

0.574780 

0.574628 

0.574627 

0.576149 

1.8 

0.597759 

0.597047 

0.596864 

0.596864 

0  598161 

2.0 

0.616944 

0.616085 

0.615944 

0.615948 

0.617074 

2.2 

0.633653 

0.632697 

0.632549 

0.632551 

0.633553 

2.4 

0.648339 

0.647267 

0.647159 

0.647161 

0.648076 

TNE,  for  both  element  models.  The  errors  are  evaluated  at  point  x  =  1.0  and  time  /  =  1.0. 
Figure  4  clearly  shows  that  convergence  cannot  be  achieved  by  increasing  the  number  of  ele¬ 
ments  alone  or  by  decreasing  the  step  size  alone. 


Convergence  For  Constant  A// Ax2 

In  the  previous  analyses  of  the  time  integration  technique  and  the  finite  element  method, 
the  results  demonstrated  that  an  increase  of  the  time  step  size  or  the  number  of  elements  alone 
does  not  guarantee  uniform  convergence.  This  phenomenon  is  similar  to  the  numerical  insta¬ 
bility  of  the  direct  finite  difference  analysis  (e.g.  [10]). 

A  modulus  M  is  proposed  here,  defined  as 


M 


A  t 

Ax2  ’ 


(52) 
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where  A/  is  the  dimensionless  time  step  size  and  Ax  the  dimensionless  element  length.  This 
modulus  M  is  of  similar  form  to  the  stability  parameter  of  the  finites  difference  method.  By 
maintaining  M  constant  for  different  values  of  TNE,  the  corresponding  values  of  Xt  are  calcu¬ 
lated  from  the  above  definition  as  follows 


Relation  between  Xx  and  Arfor  M  —  0.2 


TNE 

5 

10 

20 

30 

Results  for  constant  M 

(CE). 


NE 

Ax'1 

Ax  2 

A 1 

1 

1.0 

1.0 

0.2 

2 

2.0 

4.0 

0.05 

4 

4.0 

6.0 

1780. 

6 

6.0 

36.0 

1./180. 

given  in  Fig.  5  for  the  (LE) 


M 

0.2 

0.2 

0.2 

0.2 

and  Fig.  6  and  Table  3  for  the 


Both  sets  of  results  represent  temperture  time  histories  for  different  values  of  TNE  at 
or  —  1.0.  As  one  may  see  from  these  numerical  results,  the  convergence  of  the  finite  element 
solution  is  uniform  and  approaches  the  exact  solution  as  the  value  of  TNE  increases. 


The  error  percentage  is  show  in  Fig.  7  for  the  particular  point  at  x  —  1.0  and  t  —  1.0  for 
both  element  models.  This  figure  clearly  demonstrates  the  uniform  convergence  of  the  numeri¬ 
cal  solution  for  both  element  models  and  the  decrease  of  numerical  error  as  TNE  increases. 


Comparing  the  results  obtained  for  the  two  parts  of  this  section  it  is  concluded  that  the 
modulus  M  is  an  appropriate  parameter  for  error  control  of  the  combined  numerical  scheme  of 
finite  element  and  finite  difference  methods. 
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Fig.  6— Temperature  distribution  at  x  —  1.0  as  a  function  of  time  for  the  (CE)  model, 
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Table  3.  Numerical  results  for  the  temperature  at  x  -  1.0  for 
the  (CE)  model  with  constant  M  —  0.2. 


Time 

TNE  -  5 
A/  -  0.2 

TNE  -  10 

A /  -  0.025 

TNE  -  20 

At  -  0.0025 

TNE  -  30 

A/  -  0.0055 

Exact 

0.2 

0.069741 

0.082475 

0.107999 

.111522 

.113846 

0.4 

0.144791 

0.248910 

0.259719 

.261983 

0.263552 

0.6 

0.297175 

0.350841 

0.358799 

.360226 

0.361310 

0.8 

0.406438 

0.421731 

0.427436 

.428403 

.429195 

1.0 

0.464971 

0.473867 

0.478166 

.478931 

.479500 

1.2 

0.502502 

0.514170 

0.517715 

.518280 

.518605 

1.4 

0.534534 

0.543491 

0.549360 

.549906 

.550090 

1.6 

0.563076 

0.573147 

0.575505 

.576005 

i 

.576144 

1.8 

0.587472 

0.595604 

0.597619 

.598086 

.598161 

2.0 

0.608067 

0.614868 

1 

0.616490 

.617113 

.617074 

2.2 

0.625727 

0.631614 

0.632925 

.633571 

1 

.633553 

2.4 

0.641165 

0.646344 

0.647125 

.648091 

.648076 

A  comparison  of  the  two  element  models  shows  the  superiority  of  the  cubic  one  over  the 
linear.  This  higher  order  element  not  only  produces  a  more  stable  solution  but  also  a  much 
more  accurate  one  than  does  the  linear  element.  Even  though  the  number  of  equations  to  be 
solved  for  (CE)  is  twice  the  number  for  (LE),  and  the  computer  time  required  for  the  solution 
is  about  two  to  one,  the  (CE)  is  preferred  due  to  better  accuracy  even  for  small  TNE.  Another 
advantage  of  (CE)  is  that  any  type  of  boundary  condition  may  be  imposed  accurately  to  the 
boundary  nodal  points  since  the  generalized  coordinates  represent  both  heat  displacements  and 
temperatures. 

Similar  error  and  convergence  criteria  can  be  derived  for  the  dimensionless  ratio  A //Ax 
when  the  convection-diffusion  equation  is  solved  ( V0,K0  0)  or  when  the  pure  convection 

28 


Fig.  7 —Convergence  for  the  numerical  solution  for  M  “  0.2,  temperature  errors  at  x  —  1.0  and 
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equation  (K„  ~  0,  F,  —  1)  is  solved  Using  the  latter  as  a  test  case  we  choose  the  same 
characteristic  length  L  as  in  previous  examples  and  the  same  constants  The  boundary  condi¬ 
tion  is  assumed  to  be  a  half-sine  wave  with  period  2n  and  the  initial  conditions  are  zero,  i  e 
0(x,  0)  —  0  This  type  of  boundary  condition  is  suitable  for  investigating  solution  convergence 
since  it  has  a  continuous  form.  Results  for  the  convergence  of  the  numerical  solution  are 
presented  in  Figs.  8-11  The  first  two  figures  (8  and  9)  show  temperature  time  histories  at  *  - 
1.0  for  the  two  models  (LE)  and  (CE)  respectively  The  time  step  size  is  kept  constant  and 
TNE  is  given  four  different  values  to  show  convergence  of  the  solution  with  respect  to  Av 
Figure  10  shows  temperature  time  histories  for  the  (CE)  model  with  a  constant  value  for  TNE 
and  four  different  values  for  A t.  It  is  apparent  from  these  figures  that  uniform  convergence  is 
not  obtained  by  changing  only  A t  or  Ax.  Numerical  experimentation  showed  that  when  the 
value  of  the  ratio  A//Ax  is  maintained  constant  uniform  convergence  of  the  solution  is 
obtained.  Results  for  this  case  are  given  in  Fig.  11  for  the  (LE)  model  with  four  values  of 
TNE  and  corresponding  values  of  A /.  Similar  results  were  obtained  for  the  (CE)  model  The 
value  for  the  ratio  Ar/Axwas  equal  to  0.05. 

Comparing  the  above  results  for  the  case  of  pure  convection  it  is  apparent  that  the  ratio 
Ar/Axis  the  appropriate  parameter  for  error  control  of  the  numerical  solution 

As  a  last  test  for  the  stability  of  the  solution  experimentations  were  performed  with 
different  values  for  the  characteristic  length  L  ranging  from  1  to  5  and  constant  TNE  Since  for 
every  different  value  of  L  the  values  of  Ax  will  change,  the  ratio  Ar/Ax  and  the  modulus  M 
were  kept  constant  by  adjusting  the  value  of  A  t  accordingly.  It  was  observed  that  changes  in 
the  values  of  L  had  no  effect  on  the  solution  and  on  the  propagation  of  the  wave.  For  small 
values  of  TNE  as  L  becomes  smaller,  for  example  TNE  -  5  and  L  -  1 ,  the  numerical  solution 
becomes  more  accurate.  This  is  expected  since  Ax  is  five  times  smaller  for  L  —  1  than  for 
L  -  5  (i.e.  Ax  -  0.2  for  L  -  1  and  TNE  -  5  while  Ax  -  1.0  for  L  -  5  and  TNE  -  5). 
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Fig.  8  — Temperature  distribution  at  x  *■  1,0  as  a  function  of  time  for  the  (LE)  model,  \< 
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—  Temperature  distribution  at  x  —  1  0  or  a  function  of  time  for  the  (CE)  model,  A/-  0.025. 
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Fig.  11  —Temperature  distribution  at  x  =  1.0  as  a  function  of  time  for  the  (LE)  model,  Ar/Ax 
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c.  Numerical  Results  For  Boundary  Value  Problems. 


In  applying  the  derived  finite  element  formulation  to  diffusion-convection  problems,  the 
semi-infinite  space  is  approximated  by  the  characteristic  length  L  with  a  time-dependent  tem¬ 
perature  applied  on  its  boundary.  Three  different  cases  of  boundary  conditions  are  considered 
and  results  are  obtained  for  both  element  models. 

Boundary  conditions: 

Case  I  0(0. r)  =  1 .  i  >  0; 

=  1,  0  <  t  ^  /„, 

Case  II  8(0. t)  = 

=  0,  t  >  t; 
sin  (nt)  0  <  /  <  /„, 

CaseUI  61(0,/)  = 

=  0  /  >  /„  . 

The  boundary  condition  at  infinity  (x  =  L)  is  the  same  for  all  cases 

0(L,i)  =  0 

and  the  initial  condition  for  all  cases  is 

8  (x,  0)  =  0 


Numerical  solutions  of  the  governing  equation 


d0  I.  d9  „  d28 

b<  +  V°Tx'K°  dx2 


0 


are  obtained  by  solving  the  system  of  n  equations  represented  by 


(53) 


A,jQj  +  B,jQ,  =  Qr  (54) 

where  A,,  and  B,t  are  the  global  matrices,  given  in  terms  of  Eq.  (45)  for  the  (LE)  model  and 

Eq.  (46)  for  the  (CE)  model. 


Numerical  results  for  the  above  boundary  value  problems  were  obtained  for  the  charac¬ 
teristic  length  L  =  5,  divided  into  TNE  =  30  for  the  (LE)  model  and  TNE  =  20  for  the  (CE) 
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one.  The  corresponding  element  lengths  and  time  step  sizes  used  in  the  numerical  solution  are 
given  in  Table  4. 


Table  4.  Time  step  sizes  for  the  numerical  solution. 


TNE 

Ax 

A  / 

At/Ax 

At/A. 

LE 

30 

1/6 

0.0075 

0.045 

0.27 

CE 

20 

1/4 

0.0125 

0.05 

0.2 

1 _ J 

Temperature  time  histories  are  given  in  Figs.  (12-19)  for  the  point  at  x  =  1.0  and  tem¬ 
perature  distributions  as  a  function  of  x  are  given  in  Figs.  (20-27)  at  time  t  =  2.0.  Tempera¬ 
ture  time  histories  presented  in  Figs.  (12-15)  are  for  the  (LE)  model  and  in  Figs.  (16-19)  for 
the  (CE)  model.  Similarly,  temperature  distributions  for  the  (LE)  model  are  given  in  Figs. 
(20-23)  and  in  Figs.  (24-28)  for  the  (CE)  model.  In  each  figuie  results  are  given  for  pure  con¬ 
vection  (K„  =  0.0,  V0  =  1.0),  pure  diffusion  (Ku  =  1.0,  V„  =  0.0)  and  for  two  cases  of 
diffusion-convection,  (K„  =  0.1,  V„  =  1.0)  and  ( K0  =  1.0,  Va  =  1.0).  The  analytical  solution 
for  pure  convection  is  presented  by  a  solid  line  in  all  figures. 

For  the  first  case  of  boundary  conditions,  Figs.  (12),  (16),  (20)  and  (24),  the  numerical 
solution  shows  good  agreement  with  the  analytical  one.  The  oscillations  around  the  discon¬ 
tinuity  damp  out  as  the  wave  front  progresses.  The  error  can  be  controlled  by  the  TNE  used. 
A  finer  discretization  reduces  the  error  of  the  numerical  solution  around  the  discontinuity. 
This  finer  discretization  can  be  either  uniform  or  localized  around  the  dicontinuity.  Although 
the  TNE  used  for  both  models  is  rather  small,  the  results  obtained  depict  only  small  errors. 

For  the  second  case  of  boundary  conditions.  Figs.  (13),  (17),  (21)  and  (25),  i„  -  1.0  was 
used  which  corresponds  to  a  square  wave  propagating  through  the  half-space.  For  pure  convec¬ 
tion,  the  (LE)  model  propagates  the  square  wave  but  its  shape  is  distorted  due  to  numerical 
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dispersion.  On  the  other  hand,  the  (CE)  model  gives  a  much  better  approximation  of  the  wave 
but  the  oscillations  of  the  numerical  solution  around  the  discontinuity  are  still  present.  These 
oscillations  are  inherent  in  any  numerical  solutions  as  is  shown  by  a  Fourier  analysis  in  Appen¬ 
dix  B.  The  oscillations  and  the  error  can  be  minimized  by  a  finer  discretization  or  by  introduc¬ 
ing  some  artificial  diffusion  into  the  numerical  solution.  However,  such  an  artificially  intro¬ 
duced  diffusion  will  not  allow  a  realistic  evaluation  of  the  developed  finite  element  models,  and 
will  introduce  artificial  errors  when  true  diffusion  is  present.  From  the  obtained  results  it  can 
be  seen  that  there  is  no  phase  lag  between  the  exact  and  numerical  wave  forms,  and,  even  for 
the  rather  coarse  discretization  used  the  shape  is  well  approximated. 

The  last  set  of  boundary  conditions  represents  the  propagation  of  a  sine  wave  through  the 
half-space.  In  Figs.  (14),  (18),  (22)  and  (26)  results  are  given  for  /„  =  1.0  and  in  Figs  (15), 
(19),  (23)  and  (27)  the  sine  wave  is  continuously  applied  at  the  boundary.  For  the  (LE)  model 
results  show  some  small  error  but  the  shape  of  the  wave  is  not  distorted.  The  higher  order  ele¬ 
ment  model  (CE)  shows  an  excellent  agreement  with  the  exact  wave  even  for  the  small 
number  of  elements  used  (TNE  =  20).  For  these  boundary  conditions,  an  increase  in  the  TNE 
will  improve  the  results  for  the  (LE)  model  but  it  will  have  a  very  small  effect  to  the  alreadly 
very  accurate  results  of  the  (CE)  model.  For  all  the  cases  of  boundary  conditions  and  all 
choices  of  the  constants  K0  and  V0 ,  the  numerical  solutions  produced  accurate  results  and  the 
thermally  induced  waves  propagate  through  the  half-space  in  a  very  satisfactory  manner. 

SUMMARY  AND  CONCLUSIONS 

A  variational  formulation  for  the  convection-diffusion  equation  has  been  presented  in  this 
report,  and,  based  in  this  formulation,  two  finite  element  models  have  been  developed  for  the 
purpose  of  solving  problems  on  propagation  of  thermally-induced  induced  waves. 
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The  introduction  of  a  new  quantity,  defined  as  heat  displacement,  is  the  basis  for  a  gen¬ 
eralized  development  of  the  convection-diffusion  equation.  The  advantage  of  such  a  formula¬ 
tion  is  that  it  can  be  used  to  develop  a  displacement  formulation  for  the  finite  element  method. 
Furthermore,  due  to  the  generalized  nature  of  the  heat  displacement,  the  formulation  may  be 
extended  to  other  types  of  equations.  Another  advantage  of  the  formulation  is  the  introduced 
thermal  force,  for  which  one  should  point  out  its  significance  as  a  boundary  force. 

The  physical  conditions  for  the  semi-infinite  space  require  that  6  — 1 >  0  as  x  —  since  the 
last  nodal  point  of  the  finite  element  approximation  of  the  half  space  represents  infinity  one 
should  impose  the  above  condition  at  this  point.  The  thermal  force  is  then  zero  due  to  zero 
temperature.  This  assumption  is  not  the  correct  one  since  the  temperature  at  the  last  nodal 
point  changes  as  the  thermal  wave  propagates.  If  we  consider  the  last  nodal  point  as  a  boundary 
point  and  the  thermal  force  as  a  boundary  force,  which  is  equal  to  the  temperature  at  that 
point,  then  the  conditions  at  the  boundary  point  are  properly  adjusted.  The  presence  of  this 
boundary  force  into  the  formulation  produces  a  much  more  accurate  temperature  distribution 
close  to  the  boundary,  since  it  represents  the  effect  of  the  neglected  portion  of  the  medium. 

For  the  solution  of  the  matrix  differential  equation,  a  third  order  backward  finite 
difference  scheme  was  used  which  is  proved  to  be  unconditionally  stable.  Further,  the  conver¬ 
gence  of  the  two  element  models  was  investigated  and  some  convergence  criteria  were  dis¬ 
cussed. 

The  two  finite  element  models  developed  in  this  study  were  used  successfully  to  solve 
problems  involving  both  convection  and  diffusion  with  prescribed  boundary  conditions  for  the 
temperature.  Comparison  of  present  results  with  analytical  solutions  show  the  performance  of 
the  cubic  element  model  to  be  superior  to  the  linear  one.  However,  the  performance  of  the 
(LE)  model  should  not  be  underestimated,  especially  when  one  considers  the  crude  approxima- 
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tion  and  the  coarse  space  discretization  involved.  The  choice  between  the  two  models  for 
specific  applications  should  depend  on  the  particular  needs  of  the  problem  under  consideration 

In  conclusion,  the  generalized  form  of  the  derived  formulation,  its  efficiency  in  handling 
various  types  of  boundary  conditions,  and  its  efficiency  in  solving  not  only  diffusion-dominated 
flows  but  also  convection-dominated  flows  deserve  special  emphasis.  The  superiority  of  the 
cubic  element  model  over  the  linear  one,  especially  for  simulating  sharp  wave  fronts,  is  also 
noted.  An  extension  of  the  present  formulation  to  two  dimensional  problems  and  its  applica¬ 
tion  to  other  types  of  equations  are  planned  as  future  work. 

REFERENCES 

1.  Price,  H.S.,  Cavendish,  J.C.,  and  Varga,  R.S.,  "Numerical  methods  of  higher-order  accu¬ 
racy  for  diffusion-convection  equations,”  J.  Soc.  Petrol.  Eng.  (Sept.  1968)  pp  293-303 

2.  Lee,  R.,  Gresho,  P.,  and  Sami,  R.,  "A  comparative  study  of  certain  finite  element  and 
finite  difference  methods  for  advection  diffusion,"  1975  Summer  Computer  Sin.it'stion 
Conference,  Washington,  DC  (July  1976). 

3.  Ehlig  C.,  "Comparison  of  numerical  methods  for  solution  of  the  diffusion-convection 
equation  in  one  and  two  dimensions,"  Proc.  of  the  Int.  Conf.  on  Finite  Elements  in  Water 
Resources,  2nd,  Pentech  Press,  London,  England  (1978). 

4.  Van  Genuchten,  M.T.,  "On  the  accuracy  and  efficiency  of  several  numerical  schemes  for 
solving  the  convective-dispersive  equation,"  Pro.  of  the  Int.  Conf.  on  Finite  Elements  in 
Water  Resources,  2nd,  Pentech  Press,  London,  England  (1978). 

5.  Smith,  I.M.,  Farraday,  R.V.,  and  O’Connor,  B.A.,  "Raleigh-Ritz  and  Gaterkin  finite  ele¬ 
ments  for  diffusion  convection  problems,"  Water  Resources  Res.,  9,  3  (1973). 


55 


G  KERAMIDAS 


6.  Keramidas,  G.A.  and  Ting,  E.C.,  “A  variational  formulation  for  thermal  stress  analysis. 
Part  I:  variational  formulation,"  Nuclear  Eng.  and  Design  39,  (1976),  pp.  267-275. 

7.  Keramidas,  G.A.  and  Ting,  E.C.,  "A  variational  formulation  for  thermal  stress  analysis 
Part  II:  Finite  element  formulation,"  Nuclear  Eng.  and  Design,  39  (1976),  pp.  277-287. 

8.  Chung,  T.J.,  Finite  Element  Analysis  in  Fluid  Dynamics ,  McGraw  Hill,  Inc.,  (1978). 

9.  Salvadori,  M.G.  and  Baron,  M.L.,  Numerical  Methods  in  Engineering,  Englewood  Cliffs, 
N.J.,  Prentics-Hall  (1961). 

10.  Carnahan,  B„  Luther,  H.A.,  and  Wilke,  J.O.,  Applied  Numerical  Methods,  Clarendon  Press, 
Oxford,  2nd  edition  (1959). 

11.  Carslaw,  H.S.,  Fourier  Series  and  Integrals,  Dover  Publ.,  Inc.,  NY  (1930). 


56 


APPENDIX  A 


The  variational  formulation  in  this  study  and  the  application  of  the  finite  element  method 
produced  a  system  of  simultaneous  ordinary  differential  equations  which  may  be  represented  by 

A, 9,  4-  C„q,  -  Qr  (A.l) 

For  the  solution  of  (A.l),  the  first  derivative  is  approximated  by  a  third-order  Newton  back¬ 
ward  difference  scheme  and  at  time  step  n  has  the  form 

-  ~r  (11  q}*'  -  IS*/"-1’  +  -  2 (A. 2) 

bat 

In  terms  of  (A. 2)  the  system  (A.l)  yields 


HID,;  +  6Af Cu)qjM  -  6A tQ,M  -  -  2g/*-J,l 


[1 18, y  +  6btD«lCkJ]qfn)  -  6MD~'Qj  +  5(>tl8 q}"'"  -  9 q/n'2)  (A .3) 

+  2q/"~3)]. 

Let  the  errors  associated  with  each  time  step  be  given  by 

£/"-*\  k  -0, 1,2,3. 

If  these  errors  are  added  to  (A. 3),  one  obtains 

A0lqjM  +  EjM]  -  6AtD~tQJ(n)  +  5</Cl8(«/"-,)  +  Ef-") 

-  9 (q}"-2)  +  £/"-»)  +  2 (*/•-»  +  E}”-})))  (A. 4) 

where 

Au  -  1  iaf7  +  6AfA*'C*,.  (A. 5) 

Subtraction  of  (A.4)  from  (A.3)  yields 

AuEjM  -  8yll8£'/"-|)  -  9 £/— “  +  £)("-}>] 


£/»)  -  A,7lllS£/"~n  -  9£/"~J>  +  2£/"-})). 
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For  solution  stability  the  errors  at  various  time  steps  should  be  related  as  follows 

Efn)  <  <  £/"-2>  <  gJkE^-})  j 

with  the  amplification  matrix  g„  given  by 

g„  “  AI'B*  (A. 6)  j 

where 

I 

Bkl  -  8*J18W'°)  -  9W(U  +  2]  (A. 7)  ^ 

and 


W (2)  «£  sj  1. 

If  one  assumes  that  fV(2)  —  W(,)  —  1,  then  (A. 7)  becomes 


and  (A.6)  yields 


Bk,-n  sk/ 


8<i  “  1  i  A,j 

for  stable  solution  gn  <  8„,  then  Eq.  (A. 8)  yields 


8,  +  ±  MD,;'Ckj 


^  8, 


For  the  last  relation  to  be  valid  the  following  should  hold 


i 


(A. 8) 


I 


8,+  D~kxCkj 


8  ,j  >  1 


(A. 9) 


which  is  true  for  all  values  of  Ar  as  long  as  the  product  Dik  Ckj  is  positive  definite.  Thus,  the 


numerical  integration  scheme  is  unconditionally  stable. 


4 


As  an  example,  consider  the  case  of  the  linear  element  model  to  approximate  the  solution 
given  by  (A.l).  Assuming  that  the  characteristic  length  L ,  is  divided  in  three  equal  elements, 
the  matrices  D,j  and  C„  are  given  by 


n  _  A£ 

L/tl  x 


2  10  0 
14  10 
0  14  1 
0  0  12 
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For  a  two-point  boundary  value  problem  with 


l 

-l 

0 

0 

-1 

1 

0 

0 

K 

-l 

2 

-1 

0 

^  V, 

-1 

0 

1 

0 

0 

-1 

2 

-1 

2 

0 

-1 

0 

1 

0 

0 

-1 

1 

0 

0 

-1 

1 

the  above  matrices  reduce  to 


q(o,t)  -  q(L,t )  -  0 


-  T 


4  1 
1  4 


2K0  _ 

Ax  Ax  +  2 
Ko  Vq  2K0 

Ax  2  Ax 


Substituting  (A.  10)  into  (A.9)  one  obtains 


15Ax 


*Kq  K 
Ax  2 


Ml  ,  M 

Ax  2 


6  K0 
Ax 


K 

2 


Mi 

Ax 


JjL 

2 


>  0 


or 


(A. 10) 


which  is  true  for  all  values  of  Arand  Ax.  Therefore,  the  numerical  scheme  is  unconditionally 
stable. 
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The  presence  of  oscillation  in  a  numerical  solution  for  the  propagation  of  discontinuities  is 
known  as  Gibbs  phenomenon.  To  better  understand  this  pheomenon,  assume  that  a  rectangu¬ 
lar  wave  Hit)  of  perod  2rr  (Figure  B.l)  is  propagating  with  a  certain  velocity. 


Figure  B.l 


An  approximation  of  this  wave  can  be  obtained  by  a  Fourier  series  for  which  the  partial 
sum  of  the  first  2  n  terms  is  given  by 

Hln(t)  -  \  ±  t  TTT  «n(2*-l)/  (B.l) 

2  ir  2k—  1 

with  the  cosine  terms  all  zero.  This  partial  sum  H2„  overshoots  the  function  Hit),  as  Gibbs 
pointed  out,  by  the  amount 


1.0895  as  n  — *  oo 


In  fact,  not  only  does  this  overshoot  of  H2„  exist.  Figure  B.2,  but  the  sum  also  oscillates  about 


Hit)  with  these  oscillations  decreasing  only  away  from  the  discontinuity. 
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The  phenomenon  can  be  explained  by  rewriting  Eq.  (B.l)  as  follows 

H2At)  —  —  +  —  T  f  cos(2fc-l)x  dx 

2  n  ft  Jo 

~  +  —  f  T  cos(2A.-/)x  dx 

2  n  Jo  *- 


*.  - 1 


J_  ,  J_  C'  sin 

7  5  Jo  ti 


sin2/ix 


dx 


2  2  Jo  sinx 

The  maxima  and  minima  of  H2„  are  obtained  from  Eq.  (B.2)  by  requiring  that 

dH2nU) 


dt 


0,  7T  ^  ^  27T. 


This  requirement  is  satisfied  when 


which  is  true  for 


sin2/if 

sin/ 


0, 


(B.2) 


I  -  m  -  1,2 . 2//-1. 

2/i 

The  maxima  and  minima  alternate  and  their  values  have  been  calculated  by  Carslow  [11],  It  is 
apparert  that  the  oscillations  of  the  approximation  of  the  rectangular  wave  can  be  only  reduced 
by  including  additional  terms  into  the  partial  sum  Damping  of  the  oscillations  and 

reduction  of  the  overshoot  can  be  achieved  by  introducing  certain  artificial  factors  into  the 
terms  of  the  summation. 
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C 

PROGRAM  WAVE ( N2 . NK ) 

C 

C* 

C*************************** *********************************************** 
C* 

C*  PROGRAM  CONSTANTS 

C* 

C*  N  =  NUMBER  OF  ELEMENTS  =  NUMEL 

C*  N+l  =  NUMBER  OF  NODES  =  NNODE 

C*  WO  =  INVERSE  OF  ELEMENT  LENGTH 
C*  VO  =  DIMENSIONLESS  VELOCITY 

C*  TK  =  DIMENSIONLESS  DIFFUSIVITY 

C*  DT  =  DIMENSIONLESS  TIME  STEP  SIZE 

C*  TO  =  DIMENSIONLESS  TIME  CONST.  FOR  BOUNDARY  CONDITIONS 

C*  TMAX  =  MAXIMUM  TIME  LIMIT  FOR  INTEGRATION 
C«  I FREQ  =  CONSTANT  FOR  PRINTING  RESULTS 
C*  LCASE  =  ORDER  OF  F.E.M. 

C*  =  1  .FIRST  ORDER  ELEMENT 

C*  =  2  ,  CUBIC  HERMIT  I AN  ELEMENT 

C*  NB  =  NUMBER  OF  BOUNDARY  CONDITIONS 
C*  IBC  =  INITIAL  CONDI  TINS  CONST. 

C*  =  0  .  I.C.  SET  TO  ZERO 

C*  =  I  ,  I.C.  SPECIFIED  AT  NODES 

C*  AN  =  INTEGRATION  CONSTANTS,  N  =  0,1, 2, 3, 4, 5 
C*  R1  =  BOUNDARY  CONDITION  CONST. 

C*  RN  =  BOUNDARY  CONDITION  CONST. 

C*  Rl  =  1.0,  RN  =  0.0  SQARE  WAVE 

C*  Rl  =  0.0,  RN  =  1.0  SINE  WAVE 

C#  AT  =  GLOBAL  MASS  MATRIX 

C*  BT  =  GLOBAL  STIFFNESS  MATRIX 

C*  TE  =  TEMPERATURE  VECTOR 

C*  H  =  DISPLACEMENT  VECTOR 

C*  C  =  PARTITIONED  GLOBAL  MATRIX 

C* 

C#»**##***»*******  **##******■*##******#*•*#**#*■*#**#***#*##*#**#**«"»  **#****«# 
r# 
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COMMON/ BLOCK  1 /  N, WO, VO, TK, TO, PT, TMAX . I FREQ, LCASE, NB, IBC 
COMMON/ BLOC K2/  X 1 (200) , X2(200) ■ X3(200) , X4 (200 ) , YX ( 200 ) 
COMMON/B!  0C  K3/  X ( 200 >  ,  W  <  200 )  >  TE  ( 200 ) >  H ( 200 ) ,  FO < 200 ) 
COMMON/ BL0CK4/  NNTS, NNQS, NTS< 100) , NOS (200) 

COMMON/ BLOCKS/  AO, A1 , A2>  A3, A4 , AS, N1 , R1 , RN 
DIMENSION  AT ( N2 , N2 ) ,BT(N2,N2) , C<N1  , NK ) ,  I p ( Nf  ) 

D I MENS ION  T I ME ( 200 ) , TT ( 200 > 

D I MENS I  ON  I D 1 ( 20 ) ,  I D2 ( 20 ) 

DATA  STOP  /STOP'/ 

C 

CALL  R*STOP 

C 

C . READ  AND  PRINT  TITLE . 

r 

100  READ  (5,  120,  END=<?999>  ID1 
READ  120,  ID2 
PRINT  12S,  1 01 , 1 D2 

C 

C . READ  AND  PRINT  DATA . 

C 

READ  500,  LEASE, Rl.RN 
IF  <  LEASE.  EC1.  1  )  GO  TO  5 
WRITE ( 6 , 2015 ) 

00  TO  6 

5  WRITE (6. 2020) 

6  CONT  I NIJE 

READ ( 5 , SOS )  N.WO.VO.TK, TO, DT, TMAX, I FREQ, IBC.NB 
WRITE (6, 600)  N, WO, VO, TK, TO, DT , TMAX , I  FREQ, IBC 

ELEMENT  CONSTANTS 


NOMEL 

= 

N 

NNODE 

= 

N  +  1 

NO 

= 

LCASE* ( N  -  1) 

N1 

- 

LCA3E*NUMEL  +  1 

N2 

= 

LCASE* (NUMEL  +  1) 

MK 

= 

NK  -  1 

T 

= 

0.0 

IT 

= 

0 

. INTEGRATION  CONSTANTS  FOR  3RD  ORDER  BACKWARD  DIFF. . . . 

AO  =  11. /fc. 

A 1  =  3.0 
A2  =  ~  1 . 5 
A3  =  1./3. 

A4  =  0.0 

ICUNTE  =  0 

. READ  INITIAL  CONDITIONS 

CALL  I NCON ( N2 , NK , NUMEL , NNODE ) 

C 


*-  ' 


— TT 
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SET  BOUNDARY  CONDITIONS  FOR  TEMPERATURE . 

READ  NODE  NUMBERS  AND  SPECIFIED  TEMPERATURES 

DO  110  1  =  1,  NNODE 
I  NX  =  I 

READ (5,1 0 l 5 )  WORD , NT , TNT 
IF ( WORD. EQ. STOP)  00  TO  115 
NTS ( I )  =  NT 
TE  (  NTS  <  I  )  )  =  TNT 


110  CONTINUE 

C 

C . COUNT  NODES  WITH  SPECIFIED  TEMPERATURES . 

C 

115  NNTS  =  I  NX  -  1 

NNQS  =  N2  -  NN i S 
NX  =  NNTS 
NY  =  NNQS 


. SPECIFY  THERMAL  FORCES  AT  NODES . 

DO  130  1=1, NY 

NQS ( I )  =  I  +  1 
TE ( NQS ( I ) )  =  0.0 
30  CONTINUE 

WRITE (6, 1060) 

DO  160  1=1, NNTS 

WR I TE ( 6 , 1 065 )  I , NTS ( I ) , TE  <  NTS ( I > ) 

60  CONTINUE 

DO  15  1=1, N2 

do  is  j=  i , n:_ 

AT(I,.J)  =  0.0 
5  BT  < I , U )  =  0.0 

. FORM  THE  MATRIX  COEFFICIENT  BT(I,J> 

FROM  THE  MASS  S<  STIFFNESS  GLOBAL  MATRICES . 

CALL  ELMAT ( N2 , NK , NUMEL , AT , BT ) 


. PARTITION  THE  SINGULAR  MATRIX  BT(I,J> 

ACCORDING  TO  BOUNDARY  CONDITIONS  TO  C(I.J) 

DO  50  1=1, NNQS 
DO  50  J=1,NNQS 

C ( I , J )  =  BT ( NQS  < I ) , NQS ( J ) ) 

DO  45  K= 1 , NNTS 
5  C  <  I ,  J )  =  C  ( I ,  J )  - 

1  BT ( NQS ( I ) . NTS ( K ) ) /BT ( NTS  <  K ) , NTS ( K ) ; *BT ( NTS  <  K  > , NQS ( J ) ) 

0  CONTINUE 

CALL  DECOMP < NY, NY, C, IP, I ER) 
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C . PRINT  INITIAL  CONDITIONS  ?•<  AFTER  FIRST  TIME  STEP 


C 

c 

10 

c 

c. 

c 

20 

c 

c. 

c 


25 

C 

c. 

c 

c 

c 

c: 

c 

c 


c 


30 

C 

c. 

c 


c 


THE  RESULTS  AT  SPECIFIED  TIME  INTERVALS . 

CALL  OUTPUT ( N2 , NK , T , NUMEL . NNODE ) 

START  TIME  INTEGRATION 
T  =  T  +  DT 

PLOT  RESULTS  WITH  RESPECT  TO  TIME . 

IT  =  IT  +  1 
I WO  =  INTI  WO) 

N4  =  IWO  +  1 

I F ( MOD ( I T , I FREQ ) . NE . 0 )  GO  TO  25 
ITR  =  IT/ 1  FREQ 
TT(ITR)  =  TE(N4) 

TIME! ITR)  =  T 
CONTINUE 

FORM  THE  RIGHT  HAND  SIDE  VECTOR  CHU)  OF  THE  MATRIX  EQUATION 

C ( I , J)*X ( I )  =  Q(J> 


AND  SOLVE  FOR  X (I) . 

CALL  QVECT<N2.NX,NY,T,NUM£L,AT,BT,C, IP) 

IF  T  EXCEEDS  TMAX,  TERMINATE  I NDEGRAT I ON  ...  . 

IF(T.LE.TMAX)  GO  TO  30 
CALL  ONPLOT ( TIME. TT. ITR) 

GO  TO  100 

ICUNTE  =  ICUNTE  +  1 

PRINT  RESULTS  OR  CONTINUE  THE  SOLUTION . 

IF ( ICUNTE. NE. I FREQ)  GO  TO  20 
ICUNTE  =  0 
GO  TO  10 


65 


»ir4V- 


G.  KERAMIDAS 


C 

C . FORMAT  STATEMENTS . 

C 

1 20  FORMAT  <  20A4 ) 

1 25  FORMAT  ( 1H1 , 1 4X ,  44H***** *************************************** , 
150H* ************************** ***********************, /  15X  , 
254H*************************  ****  ******************  ******* , 

340H* ******** ******************** ***********, / 1 5X . tH****** • 

482X  >  tH*****» ,  /  15X ,  tH******  . 

520  X ,  42H  DIFFUSION-CONVECTION  IN  ONE  DIMENSION  , 20X . 
ttH****** . /15X.6H******, 

76X.50H - , 

320H - ,  tX  .  tH******,  / 15X , tH** *«**  ,  S2X  , 

9tH****** i /  15X  , 6H****»* i 1 X . 20A4» IX . tH******, /15X,tH******, 

1 1X.20A4, IX, tH******, /15X.6H******, 

282X , tH*****» , / 1 5X , 34H****** ********************* ******* . 
355H******************************************************* , 
45H*****> / 15X , 44H*********************************** ******** * , 
550H************************************************** , ////  ) 

500  FORMAT (110, 2F10. 4 ) 

505  FORMAT (15, Ft. 2. 2F3 . 4, 3F6. 4, 315) 

510  FORMAT <  6F10. 4 ) 

600  FORMAT (  •'  1 '  ,2X,  N=  '  ,  1 2 , 5X  ,  ' W0=  ' » F5 . 2  >  5X ,  ' V0=  ' ,  F8.  6, 5X  ,  'TK=  ,F8.t,5X 
1 ,  'TO*  '.F4.2.5X,  DT='  ,Ft.4,5X,  '  TMAX=' ,  Ft.  3,  5X ,  '  IFREQ= '  ,  1 3 , 5X  ,  '  B.  CON 
2D. = ' , 12.///) 

2015  FORMAT ( 40X ,  'CUBIC-HERMI T I  AN  FINITE  ELEMENT  APPROXIMATION',//) 

2020  FORMAT ( 40X , 'LINEAR  FINITE  ELEMENT  APPROXIMATION',//) 

1015  FORMAT <tX,A4, IlO.FlO.t) 

lOtO  FORMAT< 1  OX,  'BOUNDARY  CONDITIONS  FOR  TEMPERATLIE ' , 

1//.8X,  '  I  ' , 4  X ,  'NODE  ' , 5X ,  'TEMPERATURE',//) 

1 0t5  FORMAT ( 2X , 2 ( 4X , 13) , 3X , F 1 2.  3) 

9999  STOP 
END 
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SUBROUT I N£  CUBEL <  N2 , EL , NUMEL , AT , BT  > 

C** 

C******************************  ##**#*#*#****#**###*******#**#**##*#  ##*##*# 

c** 

C**  CUBEL  CALCULATES  THE  ELEMENT  MART ICES  FOR  THE 

C**  CUBIC  HERMITIAN  ELEMENT  MODEL 

C** 

C**# ****************** ***■*****#*###*###*#*#******* ****************  ********* 

c*# 

COMMON /BLOCK 1 /  N ,  WO . VO, TK , TO. DT , TMAX . I FREQ . LCASE . NB . I BC 
COMMON /BLOCK 3/  X  <  200 ) , U < 200  > , TE ( 200 > . H <  200  > , FQ 1 200 > 

DIMENSION  A<4,4) ,  B<4, A) . AT ( N2 . N2 > . BT ( N2, N2) 

C 

TKO  =  TK 
C 

DO  SO  L=l, NUMEL 

KH  =  LCASE* <L-1) 

ELEMENT  MASS  MATRIX  Ad.J) . 


A  ( 1  .  1  ) 

35 

4. *W(L>**2 

A  ( 1 . 2 ) 

= 

22. *W<L ) 

A  ( l .  3 ) 

as 

-3.*W(L)**2 

A< 1,4) 

s= 

13. *W<L ) 

A<2. 1 ) 

= 

22. *W(L) 

A<2,2) 

= 

156.0 

A<2,3> 

= 

-13. *W ( L ) 

A ( 2. 4 ) 

as 

54.0 

A (3. 1 ) 

= 

-3. *W(L)**2 

A<  3, 2  > 

S3 

- 1 3. *W  <  L ) 

A<  3, 3) 

as 

4. *W(L  >  **2 

A  <  3 , 4  ) 

as 

-22. *W  <  L ) 

A<4. 1 > 

as 

13. *W(L) 

A  (  4, 2 ) 

as 

54.0 

A  ( 4 , 3 ) 

= 

-22. *W<L) 

A<  4 , 4 ) 

156.0 

ELEMENT  STIFFNESS  MATRIX  B<I,J) . 

CONVECTIVE  TERMS . 


VI 1 

s= 

0.0 

V12 

= 

-42. #V0 

VI 3 

= 

-7. *VO*W(L) 

VI 4 

a= 

42.* VO 

V21 

=5 

42. *V0 

V22 

as 

-210.*V0/W<L> 

V23 

- 

-42. *VO 

V24 

as 

210. *V02WIL) 

V31 

as 

7. *VO*W(L ) 

V32 

= 

42.  *VC> 

V33 

= 

0.0 

V34 

* 

-42. *V0 

V41 

-42. *V0 

V42 

= 

-210. *VO/W(L) 

V43 

s 

42. *V0 

V44 

as 

210. *VO/W(L) 

-tr’ 


f 
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C 

C 


C 


c 


c 

c. 

c 

c 

c. 

c 

25 

30 


.DIFFUSION  TERMS . 

OKI  1  =  56. *TKO 

CK12  =  42. *TKO/W<  L ) 

CK13  =  -14. *TKO 
CK 1 4  =  -42. *TKO/W<L) 

CK21  =  42. #TKO/W<  L ) 

CK22  =  504.*TK0/M(L)**2 
CK23  =  42. *TKO/W(L ) 

CK24  =-504.*TK0/W(L>**2 
CK31  =  -14. #TKO 
CK32  =  42 . *TKO/W ( L  > 

CK.33  =  56. *TKO 

CK34  =  -42. *TKO/W( L ) 

CK41  =  -42. «TKO/W(L> 

CK42  =-504. #TKO/W ( L  >  **2 
CK43  =  -42. *TKO/W  <  L ) 

CK44  =  504 . *TKO/W<  L 1**2 

B<1,1)  =  VI 1  +  CK1 1 

B(l,2)  =  V12  +  CK 1 2 

B<  1 1  3)  =  V13  +  C-K 1 3 

B( 1 . 4 )  =  V14  +  CK14 

B( 2. 1 )  ■  V21  +  CK21 

B<2.2)  =  V22  +  CK22 

B (2. 3)  =  V23  +  CK23 

B<2,4)  =  V24  +  CK24 

B(3, 1 )  =  V31  +  CK31 

B(3,2)  =  V32  +  CK32 

B  <  3. 3 )  =  V33  +  CK33 

B(3,4)  =  V34  +  CK34 

B<4. 1)  =  V41  +  CK41 

B( 4 >  2 )  =  V42  +  CK42 

B( 4. 3 )  =  V43  +  CK43 

B ( 4 . 4 )  =  V44  +  CK44 

DO  25  1=1, KL 

DO  25  J= 1 , KL 

FORMATION  OF  THE  GLOBAL  MASS  MATRIX  AT(I,J) . 

AT( ( I+KL1 ) , < J+KL1 ) )  =  AT< ( I+KL1 ) , ( J+KL1 ) )  +  A(I,J> 

FORMATION  OF  THE  GLOBAL  STIFFNESS  MATRIX  BT(I,J> . 

BT<  (I-MCLl  ),  (J+KL1  )  )  =  BT <  <  I +KL 1  ) ,  ( J+KL 1  >  )  +  Bd.J) 
CONTINUE 
RETURN 
END 
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SUBROUTINE  DECOMP4NN,  NDIM,  C,  IP,  IER) 

C* 

0***»******#**********#*#***************'**##«"»*'»**'»*****«**«***«  **«**«»» 
C*  LINEAR  SYSTEM  SUBROUTINE 
C* 

C*  DECOMPOSE  THE  NN  X  NN  MATRIX  C  INTO  TRIANGULAR  L  AND  IJ  SO  THAT 
C*  L  *  U  -  F'  *  C  FOR  SOME  PERMUTATION  MATRIX  P. 

C*  MATRIX  TR I ANGULAR 1 Z  AT  I ON  BY  GAUSSIAN  ELIMINATION. 

C* 

C*  INPUT.. 

C*  NN  =  ORDER  OF  MATRIX. 

C*  NDIM  =  DECLARED  DIMENSION  OF  ARRAY  C  . 

C#  C  =  MATRIX  TO  BE  TR I ANGULAR I  ZED. 

C*  OUTPUT.. 

C*  C(I.J),  I.LE.J  =  UPPER  TRIANGULAR  FACTOR,  U  . 

C*  C(I,.J),  I.GT.J  »  MULTIPLIERS  =  LOWER  TRIANGULAR  FACTOR,  I  -  L. 

C*  I P  <  K ) ,  K.LT.  NN  =  INDEX  OF  K-TH  PIVOT  ROW. 

C*  IP ( NN )  =  (-1 )** (NUMBER  OF  INTERCHANGES)  OR  0  . 

C*  IER  =  0  IF  MATRIX  C  IS  NONSINGULAR,  OR  K  IF  FOUND  TO  BE 
C*  SINGULAR  AT  STAGE  K. 

C*  USE  SOLVE  TO  OBTAIN  SOLUTION  OF  LINEAR  SYSTEM. 

C*  DETERM ( C )  =  IP ( N ) *C < 1 , 1 ) #C < 2, 2 ) * . . . *C ( N, N ) . 

C*  IF  I P < NN ) =0 ,  A  IS  SINGULAR.  SOLVE  WILL  DIVIDE  BY  ZERO. 
C**#*******#*********##*#*-**#******#****#****************-*****  ********** 
C* 

©5  FORMAT (5X,  ■'SINGULAR  MATRIX  C  AT  STAGE  K  =',I3) 

DIMENSION  C(NDIM.NN),  IP<NN> 

C: 

C . INITIALIZE  IER,  IP<NN>  . 

C 


IER  = 

0 

IP(NN) 

=  1 

IF 

(NN  . EQ.  1)  GO  TO  70 

NM1  => 

NN  -  1 

DO 

SO  K  = 

1  ,NM1 

KP1  = 

K  +  l 

M  »  K 

DO 

10  I  » 

KP 1 , NN 

10 

IF 

( ABS ( C 

( I , K  > )  .GT.  ABS<C(M,K>)) 

I P  ( K ) 

=  M 

T  =  C (M, K ) 

IF 

(M  .EG 

.  K>  GO  TO  20 

IP(NN) 

=  -IP(NN) 

C  <  M » K ) 

=  C<  K, K ) 

C  ( K ,  K ) 

=  T 

20 

IF 

<  T  .EG 

.  0.0)  GO  TO  SO 

T  =  1. 

0/T 

DO 

30  I  = 

KP 1 , NN 

30 

C  ( I ,  K ) 

=  -C ( I , K ) *T 

DO 

50  J  = 

KP1 , NN 

T  -  C< 

M,  J) 

C(M, J) 

«  C  ( K ,  J ) 

C  <  K ,  J ) 

=  T 

IF 

(T  .EG 

.  0.0)  GO  TO  50 

DO 

40  I  = 

KP 1 , NN 

40 

C(  I > J) 

=  C(I,J)  +  C  < I » K ) *T 

30 

CONTINUE 
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60  CONT I NUE 
70  K  =  NN 

IF  <C(NN,NN)  .EG.  0.0)  GO  TO  30 
RETURN 

SO  IER  =  K 

I P  <  NN )  =  0 
WR I TE ( 6 . 85 )  IER 
RETURN 
END 


C 

SUBROUT I NE  ELMAT  <  N2 , NK , NUMEL , AT , BT ) 

C#* 

c************************************  #***■##*■»***  *********  *««««*«**  »»««-*«* 


c** 

C*#  ELMAT  FORMS  THE  OVERALL  GLOBAL  MATRIX  BT(I.J) 

C*#  IF  LCASE  =  1,  BT  < I » J )  IS  FOR  THE  LINEAR  MODEL. 

C**  IF  LCASE  =  2,  BT(l.J)  IS  FOR  THE  CUBIC  MODEL 

C#* 


C**#**#************#**####*#*##***#*  *****■**«■*****#•»*#**# **«■»**«■***#**#*** 

c** 

COMMON/BLOCK 1/  N, WO, VO. TK, TO, DT , TMAX , I FREQ. LCASE, NB, IBC 
C0MM0N/BL0CK5/  AO , A 1 , A2 , A3 . A4 , AS . N 1 , R 1 , RN 
DIMENSION  AT  <  N2 , N2 ) , BT ( N2 , N2  > 


C 


KL  =  2*LCASE 


IF (LCASE. EG. 1 )  CALL  LINEL ( N2 . KL . NUMEL, AT , BT ) 
IF (LCASE. EG. 2)  CALL  CUBEL ( N2, KL , NUMEL . AT , BT > 


DO  10  1=1, N2 
DO  10  J= 1 , N2 

10  BT(I,J)  =  A0*AT(I,J>  +  DT  *BT ( I , J ) 

RETURN 
END 
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SUBROUT  I NE  I  NIC  ON  <  N2 ,  NK  ,  NUMEL ,  NNODE  > 

C#* 

C**##*-*#*##*##**#****####*#*###*##**##***#**###**###****#******#*-nnnn,##« 
C  ## 

C**  INCON  SETS  THE  INITIAL  CONDITIONS  FOR  THE  SYSTEM 

C#*  VARIABLES  EQUAL  ZERO  IF  NB  =  0,  OR  TO  SPECIFIED 

C:**  VALUES  IF  NB  =  1. 

C*» 

c ************************************************************************* 
c#* 

COMMON /BL0CK1/  N, WO. VO, TK . TO, DT, TMAX , I FREQ . LCASE . NB . IBC 
COMMON / BLOCK 2 /  X 1 ( 200 ) , X  2  <  200 ) . X  3 ( 200 ) , X  4 ( 200 ) , Y  X ( 200 ) 

COMMON /BLOCKS/  X ( 200 ) . W ( 200 ) . TE ( 200 ) , H ( 200 ) »FG< 200 ) 

COMMON /BLOCKS/  AO, A1 , A2 . A3, A4 , AS , N1 , R1 , RN 
C 

PI  =  3.14159/2.0 
C 

DO  10  1*1, N2 

X ( I )  =  0.0 
YXC I )  «  0.0 
XI  (  I  )  *  0.0 
X2 < I )  *  0.0 
X3  < I )  *  0.0 


10  X4(I)  =  0.0 

DO  15  1=1, NNODE 
H <  I )  =  0.0 
15  CONTINUE 

C 

C: . CALCULATE  THE  LENGTH  FOR  EACH  ELEMENT 

C 

C 

DO  20  1=1, NUMEL 
20  W<I)  =  l./WO 

IF (NB. EQ. 0>  00  TO  100 

C 

C . IF  NB  *  1  SET  INITIAL  CONDITIONS . 


LO  =  N/ 10  +  1 
LI  =  N/5  +  1 
H<L0-1)  =  -W<LO) 

DO  25  I=L0.L1 
TE< I )  =  1.0 

25  H ( I )  =  H< 1-1 )  +  W( I ) 

H(LO-l)  =  0.0 
DO  30  I=L1 » N 
30  HI  1  +  1 )  =  H(L1 > 

C 

IF (LCASE. EQ. 2)  GO  TO  40 
DO  35  1=1, NNODE 
35  X  < I >  =  HI  I) 

GO  TO  50 

40  DO  45  1*1, NNODE 

X  <2*1-1 )  =  TE< I ) 

45  XI2*I)  *  H ( I ) 

50  DO  55  1*1, N2 

XIII)  *  X  <  I  > 

X2< I >  =  X 1  ( I ) 

X3( I )  =  X2< I ) 

55  X4  < I )  =  X3 ( I ) 

C 

100  RETURN 
END 

71 


r 


►*►*000000  o  o  o  o  o  o  o  o  o  o 

(NO**  •  •  • 


G  KERAMIDAS 


SUBROUT  I NE  L I NEL ( N2 ,1  L. NUMEL , AT . BT ) 

C** 

c  **************************************************************************** 
c** 

C**  LI NEL  CALCULATES  THE  ELEMENT  MATRICES 

C**  FOR  THE  LINEAR  ELEMENT  MODEL 

C** 

c*********************************************************************** 

c** 

COMMON/ BLOCK 1/  N, WO, VO, TK, TO, DT, TMAX , I FRED, LCASE, NB, IBC 
C0MM0N/BL0CK3/  X  I  200 ) , W ( 200 ) , TE ( 200  > , H ( 200 ) . FG ( 200 ) 

DIMENSION  A ( 2, 2 ) . B ( 2 . 2 ) ,AT(N2,N2) ,BT(N2,N2) 


DO  1 5 

L=  1 

, NUMEL 

KL1 

= 

LCA 

SE*(L- 

ELEMENT  MASS 

MATRI 

A  (  1 

,  i  > 

= 

2.0 

A  (  1 

,2) 

= 

1 .0 

A  ( 2 

,  1  ) 

= 

1.0 

A  <  2 

,2) 

= 

2.0 

- ELEMENT  STIFFNESS  MATRIX  B<I,J> 


TKO  =  TK 
WO  =  l./W<L> 

_ CONVECTIVE  TERMS . 

Vll  =  -3.*W0*V0 
V12  =  3.*WO*VO 
V21  =  -3. *WO*VO 
V22  =  3. *WO*VO 

_ DIFFUSION  TERMS . 

CKl 1  =  6. *TK0*W0**2 

CK12  =  -6. *TK0*W0**2 
CK21  =  -6. *TK0*W0**2 
CK22  =  6. #TK0*W0**2 

B  (  1 ,  1  >  =  Vll  +  CK11 
B< 1 >  2)  =  V12  +  CKl 2 
B(2, 1 )  =  V21  +  CK21 
B  <  2, 2  >  =  V22  +  CK22 
DO  10  1=1, KL 
DO  10  J=1,KL 

_ FORMATION  OF  THE  GLOBAL  MASS  MATRIX  AT(I.J) . 

AT ( ( I +KL 1 > , ( J+KL 1 ) )  =  ATI (I+KL1 ) , (J+KL1 > )  +  AII.J) 

_ FORMATION  OF  THE  GLOBAL  STIFFNESS  MATRIX  BT(I,J) . 

BTI < I+KL1 ) , ( J+KL1 ) )  =  BTI < I+KL1 > , < J+KL1 > >  +  BII.J) 
CONTINUE 
RETURN 
END 
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SUBROUT I NE  OUTPUT ( N2 , NK , T , NUMEL . NNODE > 

C** 

c*********************************************************************** 

c** 

C**  OUTPUT  PRINTS  RESULTS  AT  SPECIFIED  TIME 

C**  INTERVALS  BY  I  FREQ . 

C**  OUTPUT  CALLS  EXTERNAL  SUBROUTINE  ONPLOT (  ) 

C**  FOR  PLOTING  RESULTS. 

C** 

C*«****«**»**««**««**»«***«**«««*««*«***e*******«****««***««*«#**#*»**** 
C  ■** 

COMMON/BLOCK  1 /  N, WO. VO. TK, TO. DT, TMAX , I  FREQ . LEASE . NB . IBC 
COMMON/ BLOC K2/  XI (200) . X2<200) . X3(200> , X4(200) . YX(200) 

COMMON /BLOCK 3/  X ( 200 > . W ( 200 ) . TE ( 200 ) >  H ( 200 ) , FQ< 200 ) 

COMMON/ BLOCKS/  AO. A1 . A2. A3, A 4 , A5 > N I , R1 , RN 
DIMENSION  E'H < 200 > .  DX  (200)  ,  RN0DE(200> 

C 

C 

C . COMPUTE  THE  DISPLACEMENT  !<  TEMPERATURE  FOR  THE  LINEAR  ELEMENT . 


TKO  = 

TK 

SUM1  = 

0 . 0 

SUM2  = 

0.0 

IF (LCASE. 

EQ.2) 

GO  TO  100 

DO  10  1=  1,  NNODE 
10  H  (  I  )  =  X  (  I  ) 

DO  15  1=1, NUMEL 

15  TE < I )  =  ( H < I  +  1 )  -  H(I))/W(I) 

DO  20  1=1, NUMEL 

DH  (  I  )  =  AO*  < H ( I  )  +  H<  1  +  1  )  )  -  (  YX  (  I  )  +  YX ( I  + 1  )  ) 

DX ( I )  =  2. * ( H ( I )  -  H< 1+1 ) >*VO*DT/W< I ) 

SUM1  =  SUM 1  +  DH ( I ) 

20  SUM2  =  SUM2  +  DX (  I  ) 

SUM 3  =  ABS ( SUM 1 )  -  ABS<SUM2> 

C 

00  TO  150 
C 

C . COMPUTE  THE  DISPLACEMENT  S<  TEMPERATURE  FOR  THE  CUBIC  ELEMENT . 

C 

100  DO  110  1=1, NNODE 

TE ( I )  =  X< 2*1-1) 

110  H( I )  =  X  <  2* I ) 

DO  115  1=1, NUMEL 

DH < I )  =  <TE( 1  +  1  )  -  TE( I ) )*A0  -  YX(2*I  +  1)  +  YX(2*I-1> 

1  -  6 . * ( ( H ( I + 1 )  +  H ( I > ) *A0  -  YX (2*1+2)  -  YX (2*1 ) ) /W ( I ) 

DX ( I )  =  1 2 . *DT * ( VO* ( H ( I + 1 ) -H ( I ) )  +  TKO*(TE( 1+1 )-TE< I ) ) )/W< I )**2 
SUM  1  =  SUM1  +  DH  < I  ) 

115  SUM2  =  SUM2  +  DX ( I ) 

SUM3  =  ABS ( SUM 1)  -  ABS(SUM2) 

C 

C . PRINT  THE  RESULTS . 

C 

150  WR I TE (6,610)  T 

WR I  TE  ( 6 , 620 )  SUM  1  ,  SIJM2 ,  SUMS 

WRITE (6, 605 ) 

DH( NNODE)  =  0.0 
DX( NNODE)  =  0.0 

WRITE (6, 61 5)  < I , H ( I ) , TE ( I ) , DX  < I ) , DH ( I > , 1  =  1 .NNODE) 
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C . PLOT  RESULTS . 

C 

DO  SO  1=1, NNODE 
SO  RNODE(I)  =  FLOAT  (  I  ) 

C 

C ALL  ONPLOT ( RNODE , TE , NNODE > 

C 

600  FORMAT ( 5  X . E 1 4 . 6 ) 

60S  FORMAT  (////,  NODAL  POINT  H-Dl SPLAT EMENT  TE-TEMPERATURE 
1  DH/DX  DH/DT  , //> 

610  FORMAT  (  1  HI  .  5X  ,  TIME  =  .  F6.  3,  /  /  > 

6 1 5  FORMAT ( 1  I  1 2 , 4  X , 4F 1 5. 6 ) 

620  FORMAT ( 20 X , ' DH/DT  =  '  , F 1 0. 6 , 5X ,  DH/DX  = ' , F 1 0. 6 , 5X , F 1 2. 6 . / / > 
RETURN 
END 


SUBROUT  I NE  OVEC T (N2, NX , NY, T, NUMEL , AT , BT , C , IP) 

r  *  * 

C  **«****«--*■»***  **************************************************************** 

c** 

C**  OVECT  SETS  THE  RIOHT  HAND  SIDE  VECTOR  0 ( J > 

C**  AND  FORMS  THE  SYSTEM  OF  EQUATIONS 

C**  C(I,J)*X(I)  =  Q(.J>  ,  I... I  =  1  i  NY 

C**  QVECT  CALLS  SUBROUTINE  SOLVE  FOR  THE 

C**  SOLUTION  OF  THE  SYSTEM  OF  EOS. 

C** 

c************* ********************* ***********************«^***** ******* 
c** 

COMMON /BLOCK  1 /  N, WO, VO , TK, TO, DT , TMAX , I  FREQ, LEASE , NB , IBC 
COMMON / BLOCK 2/  X 1 ( 200 ) , X 2 ( 200 >,X3< 200 ) , X 4 ( 200 > , Y X ( 200 ) 

COMMON/ BLOCKS /  X ( 200 ) , W ( 200 ) , TE ( 200 ) , H ( 200 ) , FQ  < 200 > 

COMMON/ BL0CK4/  NNTS, NNQS , NTS ( 100) ,NQS(200) 

COMMON/ BLOCKS/  AO , A 1 , A2 , A3 , A4 , AS , N 1 , R 1 , RN 
DIMENSION  AT ( N2 , N2 ) , BT < N2 , N2 ) , C < NY , NY ) , IP(NY) 

D I MENS ION  F ( 200 ) , Q < 200 ) , L I  ST X ( 1 00 ) , L I STY ( 200 > 

C 

FI  =  3.14159 
C 

C . BOUNDARY  COND I T I ONS . 

C 

IF(T.LE.TO)  GO  TO  10 
TE ( 1 )  =  0 . 0 
GO  TO  15 

10  TE ( 1  )  =  R 1  +  RN*SIN ( T  *PI ) 

15  CONTINUE 

IF1NB.EQ. 1  )  TE (  1  )  =  0.0 


74 


SKI  Ml  Ml  IK  \  SI  H  M  Rl  PORI  4I:< 


r 

C . FORM  THE  RIGHT  HAND  SIDE  VECTOR  0(1) . 

C 

DO  20  1=1, N2 
F  Q ( I )  =  0.0 

20  YX  <  I  )  =  A1 *X 1 < I )  +  A2*  X2 (  I  >  +  AH*x;-:(I)  +  A4*a4(I> 

C 

IF  (LEASE.  EG!.  2)  GO  TO  25 

FG  (  1  )  =  -4..*TK«DT*TE<  1  >/W<  1  ) 

FG'(Nl)  =  6.  *TK*DT»TE(N)/W(N> 

GO  TO  30 
25  CONT I NUE 

FQ ( 2 )  =  -420. *TK*DT«TE< 1 ) /W( 1 > 

F0(N2)  =  420 .  *TK*DT *TE  <  N+ 1  )  /W ( N  > 

30  CONTINUE 

C 

DO  40  1=1, N2 
F ( I )  =  FQ  < I ) 

DO  40  J= 1 , N2 

40  F <  I  )  =  F (  I  )  +  AT(  I,,U*YX(.J) 

C 

DO  50  1=1, NY 

Q  < I )  =  F ( NQS ( I ) ) 

DO  45  K= 1 , NX 

45  Q  (  I  )  =0(1)  -  BT  <  NG'S  (  I  )  ,  NTS  ( K  >  >  /BT  (  NTS  ( K )  .  NTS  ( K  )  )  *F  ( NTG  ((■')> 

50  CONTINUE 

C 

C . SOLVE  THE  SYSTEM  OF  EQUATIONS . 

C 

CALL  SOLVE  ( NY ,  NY ,  C ,  C! ,  I P ) 

C 

DO  1 00  I = 1 . NY 
1 00  X ( NQS ( I ) )  =  Q ( I ) 

IF  (  LCASE.  EG!.  2 )  00  TO  120 
X  < 1 )  =  X ( 2  >  -  TE< 1 >*W( 1 ) 

I F  (  IBC..NE.  1  ) 

*  X ( N+ 1 )  =  X ( N )  +  TE ( N ) *W( N  > 

GO  TO  130 
120  CONTINUE 

DO  125  1=1, NX 

125  X  <2#NTS( I >  —  1 )  =  TE(NTS( I ) ) 

130  CONTINUE 

DO  150  1=1, N2 

X  4  <  I  )  =  X  3  (  I  ) 

X3 < I )  =  X2(  I ) 

X2( I )  =  XI ( I  ) 

X I  <  I  )  =  X  <  I  > 

150  CONTINUE 

RETURN 

END 
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c 

SUBROUTINE  SOLVE  (NN,  NO I M ,  C,  0,  IP) 

C  * 

C  t*«****»JHt«*»*#*t(*»«****t**tt******«»t*«*M**«**«*n*(t««(j»(mHHHt*  **■»**« 
C*  SOLUTION  OF  LINEAR  SYSTEM,  C*X  =  Q  . 

C*  I NPUT . . 

C*  NN  =  ORDER  OF  MATRIX. 

C*  NDIM  =  DECLARED  DIMENSION  OF  ARRAY  C  . 

C*  C  =  TRIANGULARIZED  MATRIX  OBTAINED  FROM  DECOMP. 

C*  0  =  RIGHT  HAND  SIDE  VECTOR. 

C*  IP  =  PIVOT  VECTOR  OBTAINED  FROM  DECOMP. 

C*  DO  NOT  USE  IF  DECOMP  HAS  SET  IER  .NE.  0. 

C*  OUTPUT . . 

0*  0  =  SOLUTION  VECTOR,  X  . 

C**  ********************  ********#*■»********•******■»■»******■»*•**«•***■*«■*  *■**#* 

c* 

DIMENSION  C(NDIM.NN)  ,  Q < NN ) ,  IP(NN) 

C 

IF  (NN  .EG.  1)  GO  TO  50 
NM1  =  NN  -  1 
DO  20  K  =  1 , NM1 
KP1  =  K  +  1 
M  =  I P  (  K ) 

T  =  Q  ( M  ) 

Q  ( M )  =  G  ( K ) 

Q  ( K )  =  T 

DO  10  I  *  KPl.NN 
10  0(1)  =0(1)  +  C< I > K)*T 

20  CONTINUE 

DO  40  KB  =  l.NMl 

KM1  =  NN  -  KEi 

K  =  KM1  +  1 
Q  ( K )  =  0.  ( K )  /C-  ( K ,  K ) 

T  =  -Q  ( K  > 

DO  30  I  =  1 , KM1 

30  0(1)  =  Q(I)  +  C(I,K)*T 

40  CONTINUE 

50  0(1)  =  Q< 1 )/C< 1 , 1 ) 

RETURN 

END 


